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Abstract 

Treatment with high energy ionizing radiation is one of the main methods in modern 
cancer therapy that is in clinical use. During the last decades, two main approaches to dose 
calculation were used, Monte Carlo simulations and semi-empirical models based on Fermi- 
Eyges theory. A third way to dose calculation has only recently attracted attention in the 
medical physics community. This approach is based on the deterministic kinetic equations of 
radiative transfer. Starting from these, we derive a macroscopic partial differential equation 
model for electron transport in tissue. This model involves an angular closure in the phase 
space. It is exact for the free-streaming and the isotropic regime. We solve it numerically 
by a newly developed HLLC scheme based on [8], that exactly preserves key properties of 
the analytical solution on the discrete level. Several numerical results for test cases from the 
medical physics literature are presented. 

1 Introduction 

Mathematical methods play an increasing role in medicine, especially in radiation therapy. Sev- 
eral special journal issues have been devoted to cancer modeling and treatment, cf. flU |5l [6l \Y2\ 
among others. 

Together with surgery and chemotherapy, the use of ionizing radiation is one of the main 
tools in the therapy of cancer. The aim of radiation treatment is to deposit enough energy in 
cancer cells so that they are destroyed. On the other hand, healthy tissue around the cancer cells 
should be harmed as little as possible. Furthermore, some regions at risk, like the spinal chord, 
should receive a dose below a certain threshold. 
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Most dose calculation algorithms in clinical use rely on the Fermi-Eyges theory of radiation 
which is insufficient at inhomogenities, e.g. the lung. This work, on the other hand, starts with a 
Boltzmann transport model for the radiation which accurately describes all physical interactions. 

Until recently, dose calculation using a Boltzmann transport equation has not attracted much 
attention in the medical physics community. This access is based on deterministic transport 
equations of radiative transfer. Similar to Monte Carlo simulations it relies on a rigorous model 
of the physical interactions in human tissue that can in principle be solved exactly. Monte Carlo 
simulations are widely used, but it has been argued that a grid-based Boltzmann solution should 
have the same computational complexity [9]. Electron and combined photon and electron radia- 
tion were studied in the context of inverse therapy planning, cf. ff35ll34ll and most recently 11361 . 
A consistent model of combined photon and electron radiation was developed [23 that includes 
the most important physical interactions. Furthermore, several neutral particle codes have been 
applied to the dose calculation problem, see f[T9ll for a review and most recently If38l . 

In this paper, we want to study a macroscopic approximation to the mesoscopic transport 
equation. After the problem formulation in section [2l we derive the approximation of the macro- 
scopic model in section [3l This approximation consists of a system of nonlinear hyperbolic par- 
tial differential equations, whose properties we briefly discuss. Due to the possibility of shock 
solutions, hyperbolic PDEs have to be solved with great care. In section[5l we introduce a scheme 
which is adapted to the problem at hand. Numerical results for tests from the medical physics 
literature are presented in section [6] 

2 A deterministic model for dose calculation 

A ray of high energy electrons that interacts with human tissue is subject to elastic scattering 
processes and inelastic ones. It is this latter process that leads to energy deposition in the tissue 
i.e. to absorbed dose. 

To formulate a transport equation for electrons we study their fluence in phase space. Let 
y/(r, e, Q.) cos QdAdCldedt be the number of electrons at position r -r being a vector in 2D or 3D 
space- that move in time dt through area dA into the element of solid angle dQ. around Q. with 
an energy in the interval (e,e + de). The angle between direction Q. and the outer normal of dA 
is denoted by 0. The kinetic energy e of the electrons is measured in units of m e c 2 , where m e is 
the electron mass and c is the speed of light. 
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2.1 Boltzmann transport equation 

The transport equation can generally be formulated as 03] 



Q-Vy(r,e,ft) = p in (r) J J a- m (e' ,e,Q! -a)\ir{r,E' ,Q.')d€l'de' 

+ p el (r) J a el (r,e,a'-a)xi/(r,e,a')dQ! 

y 2 

- Pin (rX(e)^(r,e^) 

- p d (r)o^(r,e)vr(r,e,Q), (1) 

with o m being the differential scattering cross section for inelastic scattering, and a e i the dif- 
ferential cross section for elastic scattering; o^° l = J S 2 o m dQ. and o^ 1 = J s i o e \dQ. are the total 
cross sections for inelastic and elastic scattering, respectively; p m and p e j are the densities of the 
respective scattering centers. 

Explicit formulas for the cross sections that we used in this model can be found in section 
12.31 They are based on the model developed in [|22l . The energy integration is performed over 
(e,°°) since the electrons lose energy in every scattering event. Also, we consider only electron 
radiation. Equation (OQ) could also be used to model electrons which are generated by the inter- 
actions of photons with matter, as in [|22l . In this case we would have an additional source term 
on the right hand side for the generated electrons. 

Besides the transport equation one needs an equation for the absorbed dose. It was derived 
in [|22l as an asymptotic limit of a model with a finite lower energy bound £ s > 0. The formula 
is exact if one chooses the lower energy limit £ s = 0, as we do here. 

p(r) Jo 

with 

xif^(r,e):=J y{r,£,Q!)dQ! , 

s 2 

T being the duration of the irradiation of the patient and p the mass density of the irradiated 
tissue. If all quantities are calculated in SI units, equation © leads to SI units J/kg or Gray (Gy) 
for the dose. 

S is the stopping power related to the inelastic cross section. It is defined as 



S(r,e) =pinM J e'o- m {e,e')de'. 



2.2 Continuous slowing-down approximation 

Electron transport in tissue has very distinctive properties. The soft collision differential scatter- 
ing cross sections have a pronounced maximum for small scattering angles and small energy loss. 
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This allows for a simplification of the scattering terms in the Boltzmann equation. The Fokker- 
Planck equation is the result of an asymptotic analysis for both small energy loss and small 
deflections. It has been rigorously derived in [31] and has been applied to the above Boltzmann 
model in [[22] . However, some electrons will also experience hard collisions with large changes 
in direction and energy losses which have to be described by Boltzmann integral terms. Thus 
we only use an asymptotic analysis to describe energy loss, called continuous slowing-down 
approximation. This approximation has a greater domain of validity than the Fokker-Planck ap- 
proximation. The Boltzmann equation in continuous slowing-down approximation (BCSD) is 

+ p el (r) J a A (r,e,£i! 'Cl)Y(r,e,a')d£i! 

y 2 

- pinMcVtoKeM^e^) 

+ ^{S(r,£) ¥ (r,£,a)) (3) 



with 



POO 

°in SD = / oin(e,e',Ai)de'. 
Jo 



A truncation in the energy space is introduced, that does not allow particles with arbitrary high 
energy, 

lim \jf(r,e,Cl) = 0. (4) 

e — >oo 

In the numerical simulations, we use a sufficiently large cutoff energy. Furthermore, we prescribe 
the ingoing radiation at the spatial boundary, 

\ff(r,e,Q.) = yb{r,e,Q.) for n-^<0, (5) 

where n is the unit outward normal vector. 



2.3 Modeling of Scattering Cross Sections 
2.3.1 Henyey-Greenstein Scattering Theory 

The detailed interactions of electrons with atoms give rise to complicated explicit formulas for 
the scattering coefficients. Because of this, many studies use the simplified Henyey-Greenstein 
scattering kernel for elastic scattering 0, 

The parameter g, which can depend on r, is the average cosine of the scattering angle and is 
a measure for the anisotropy of the scattering. The case where g < 1 matches an anisotropic 
scattering configuration. 
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2.3.2 Mott and M0ller Scattering 

A more realistic model for elastic and inelastic scattering of electrons in tissue has been devel- 
oped in II221 . This model introduces material parameters (namely densities p e and p c , ionization 
energy £g and effective atomic charge Z). The energy integration for inelastic scattering is cut-off 
at e B . 

The model uses the Mott scattering formula for elastic scattering of an electron by an ion 

mm 



~n Z 2 (r)rf(l+e) 
OM 0tt (r,e,n'-n)- 



4[e(£ + 2)] 2 (l + 277(r,£)-costf) 2 



1 £ £ + 2 . 2 $ 

1 - 77 vf Sln 7-T 

(l + £) 2 2 



with # = arccos(fy • CI), and £ is the outcoming electron energy in m e c 2 units. Here, a ~ 1/ 137 
is the fine structure constant, Z is the atomic number of the irradiated medium, r e is the classical 
electron radius. Z depends on r to account for heterogeneous media. To avoid an otherwise 
occurring singularity at # = a screening parameter 

, . % 1 a 1 Z 1 l\r) 
rj(r,e) = 



£(£ + 2) 

can be introduced [|39ll that models the screening effect of the electrons of the atomic shell, 
denoted by a. 

The inelastic scattering process is M0ller scattering, where an electron impinges an atom that 
releases itself an additional electron 

e~ + a — >• a + + 2 e~ 

For this process, the electrons can be considered indistinguishable. The electron which has the 
higher energy after the collision is called primary electron, the other electron secondary. Due 
to kinematical reasons of the scattering processes the range of solid angles in M0ller scattering 
is restricted. After the collision, the angle between the directions of the electrons is at most 
n/2. For an angle in [0, it /A], the electron with energy £ is the primary electron, for an angle in 
[7l/4, 7t/2], it is the secondary electron. Therefore the M0ller cross section can be written as 

° M = ^ M ^{0<f2-fl'<V2/2} + ^M,5Z{V2/2<f2-f2'<l}' 

where % denotes the characteristic function of a set, 

c M {£',£,Ci! CI) = a M (£ / ,£)5 M (Ai,Aip)^ ) M =Cl' Ci, 
is the M0ller differential cross section of primary electrons and 

o M CI' -CI) = a M (e\e)5 M ^^ s )^, n = Cl'Cl, 
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is the M0ller differential cross section of secondary electrons. Here, 



. , , 2nrl{e' + \) 
c7 M (e,e)- 



£'(£' + 2) 



1 



1 1 

+ ■ 



2e'+l 



and 



£ 2 ( e /_ e )2 ( e / + i)2 ( e / + i)2 e ( e /_ e ) 

(e'-e B ) 



££' + 2 

i 7 e + 2 

££' + 2 

e 7 e + 2 



for e > 



e < 



2 

(£'-£b) 



In the simulations the model parameters p e i, pi n , £b and Z are fitted to tabulated values taken 
from the database of the PENELOPE Monte Carlo code (33. 



3 Partial Differential Equation Model 

We will try to reduce the cost of solving system (OQ) by assuming a minimum entropy principle 
for the angle distribution of particles. This principle has been first proposed by Jaynes [|24l as 
a method to select the most likely state of a thermodynamical system having only incomplete 
information. It has subsequently been developed in 11281 . [|27l . 0J and lfT5ll . among others, and 
has become the main concept of rational extended thermodynamics ll30l . A full account and an 
exhaustive list of references on the historical development can be found in ETI . 
We define the first three moments in angle: 

v/°)(r,e)= [w(r,£,a)da, (7) 



yz-Wfoe^ / Q\ir(r,e,Q)dQ, (8) 
Js 2 

y/i 2 ) ( r , e) = / (Q <g> Q) y(r, £, a)dQ, (9) 
Js 2 

where we note that y/°) is a scalar, y/ 1 ) is a vector and i// 2 ) is a tensor. 

If we integrate the system © over Q., we can derive the following equations, 

V,^ = ^V (0) ), (10a) 

v,v/( 2 ) = -(r M + r M o,t)^ (1) + ^(^ (1) )- (iob) 

We have introduced the transport coefficients 

Ae-e B )/2 rl 

r in (r,e) = 7tp in (r) / / (l-n)o m (£,£ f ^)d^d£\ (11) 

Jeb J— I 

T e i(r,£) = Kp el (r) J (l-y,)o & \(£,nW- O 2 ) 
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(a) Eddington factor. (b) Eigenvalues. 

Figure 1: Eddington factor % and system eigenvalues versus anisotropy parameter \a\ 



These coefficients and the stopping power can be computed for both Henyey-Greenstein and 
Mott/M0ller scattering. Explicit expressions can be found in 11221 [FT! . 

The remaining problem is the computation of moment y/^ as a function of and l/A 1 ). 
The Minimum Entropy M\ closure for electrons ifTOl can be derived in the following way. To 
close the system we determine a distribution function that minimizes the entropy of the 
electrons, 



s 2 



V/log ydQ., 



(13) 



under the constraint that it reproduces the lower order moments, 

w M£ dQ=v (0) and / Q.w ME dQ. = w w . (14) 

Is 2 Js 2 

By using this entropy, we have implicitly assumed that the electrons obey classical Maxwell- 
Boltzmann statistics. This is justified, since here quantum effects can be neglected. 

Analogous to the calculations in [27J we can show that the entropy minimizer has the follow- 
ing form, 

¥me = floexp(-Q-ai), (15) 

where ciq is a non-negative scalar, and a\ is a three component real valued vector. This is a 
Maxwell-Boltzmann type distribution and ciq, a\ are (scaled) Lagrange multipliers enforcing the 
constraints. An important parameter is the anisotropy parameter a, 



a 



¥ 



(i) 



¥ 



whose norm is by construction less than or equal to one. If we compute the different moments 
of the distribution function given by (fT5l) we obtain, 



(o) 



sinh(|ai|) m 



fli | 



Anao 



sinh( \a\ )(1 




a i 


coth(|fli|)) 




a \ 


3 



(16) 
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In fact, these relations can be combined to give, 

a = 1 Vl [/ a h (17) 



1 - 


a\ 


coth(|ai|) 




CL\ 


2 



or by taking the modulus, 

\a\ = "vi~n/ !_ (lg) 





coth(|ai|) - 1 









The relation (fT8l) cannot be inverted explicitly by hand, i.e. we cannot express \a\ | as a func- 
tion of a in a closed form. However, this relation determines a unique solution which can in 
principle be computed. If we assume that we know a\, y/ 2 ) can be computed as 



where 

x = — — "'"W'" 1 " ' " (20) 



\a\ 


2 -2| 


a i 


coth(|ai|)+2 




Cl\ 


2 



is a function of a by means of (|T8l) . 

For its efficient numerical evaluation, the Eddington factor has to be approximated. Several 
possibilities exist: 



• One could solve the closure relation (1181) for \a\\ e.g. by a Newton iteration in each step 
during the simulation. 

• One could precompute a table that gives the Eddington factor % as a function of a. 

• One could approximate %{(X) by a suitable special function. 

The second approach has been followed in [fTTll . It is advantageous only if the space in which 
one interpolates is low-dimensional. For more moments, this approach becomes more expensive, 
and the first approach appears to be more advantageous. 

In some cases, an ansatz for % can provide a good approximation. This is the approach we are 
following here. The Eddington factor % can be approximated by a very simple rational function, 

, . a 6 a 6 + a4Ct 4 + a2Ct 2 + ao 
X{a) " + + bo. <2l) 

This approximation is very accurate (the difference with exact curve is about 1(T 15 ). The coef- 
ficients are given by 

a = 0.762066949972264, b = 2.28620084991677, 

a 2 = 0.219172080193380, Z? 2 = -2.10758208969840, 

a 4 = -0.259725400168378, ( } 

a 6 = 0.457105130221120. 
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4 Properties of the System 



In the literature, the system that has been thoroughly investigated (both analytically and nu- 
merically) is system (fTOl) restricted to its conservative terms, without external sources, but with 
time-dependence . 

In the present work, we adapt a pseudo-time technique. We focus on the spatial discretization 
and use a standard discretization for the terms on the right-hand side. Thus we consider 

— i/^+V^v/ 1 ) =0, (23a) 
at 

^)+vM 2) (w { °\w {l) )=^ (23b) 

with the closure (fT8~l) . 

The Eddington factor % is shown in Figure [U Furthermore, we show the system eigenvalues 
in two dimensions. In the isotropic regime (anisotropy parameter zero), they coincide with the PI 
eigenvalues. On the other hand, in the case of free- streaming {\a\ = 1), they coincide and have 
absolute value one. Thus the system (fTOl) is hyperbolic and the speed of propagation is limited 
by one. Moreover the system is hyperbolic symmetrisable [fT51 . 

System (1231) closed by the relation (fT9l has been analyzed thoroughly in 11231 . There, so- 
lutions to Riemann problems are constructed and invariant regions are computed. Since the 
reconstruction (fT5T) of the kinetic distribution \y is always positive, it can be expected that system 
(fT9l) - (l23l) must admit a positive solution y/°) and a limited flux ||a|| < 1. To our knowledge, 
however, there exists no proof of this fact. The invariant regions computed in ||23l only cover a 
subset of all admissible values. For a related model [18J, bounds were proved, but only in ID 
and steady state. Nevertheless, we construct a scheme which preserves exactly the positivity of 
y/°) and the flux limitation, i.e. the convex set of the admissible states of the system (1231) is (HI 

^={(y (0) y i) ): V (0) >0, |V (1) |< V (0) }- 

In the absence of sources or boundaries, the total mass, momentum and energy are conserved. 

In addition, the minimum entropy system recovers the equilibrium diffusion regime as a 
relaxation limit for large absorption coefficients [fT3l . 

In a two- or three-dimensional geometry, we have in addition [81 : Let n be the unit normal 
vector to an interface; then the system exhibits two acoustic waves, with velocities and 
Aft(n), supplemented by a contact wave with velocity /3(n). The quantity /3 ■ n satisfies the fol- 
lowing inequality Xl(h) < j8(n) -n < A^(n). The Riemann invariants associated with the contact 
wave are {/3,I1}. They are defined by the relations 

Vi = (II+yb)j3, (24a) 
V/ 2 = (^o + n)/3®/3+n/J. (24b) 
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5 Numerical Method 



The properties of the continuous model should be reproduced by the numerical scheme. In 
particular the positivity and flux limitation constraints are fundamental. An HLL scheme [|20l 
can be constructed [HI [H, that satisfies the required properties. However such an approach 
cannot capture the contact discontinuity. To prevent this failure, an HLLC scheme [3] has been 
derived, that resolves the contact discontinuity and satisfies the physical constraints. 

To complete this presentation of the numerical approximation, we mention that a suitable 
high order extension that preserves both the positivity and the flux limitation can be derived, 
relying on an appropriate limitation procedure. 



5.1 An HLL scheme for the free transport M\ angular moment system 

In this section, we derive a Finite Volume method, issued from the HLL method EDI to discretize 
the free transport equation (l23l) . Put in other words, we omit the source terms and we consider 
the one dimensional generic conservative system 

JUr + j^02O]=O, (25) 



where 

V vi 

and & stands for the flux of the M\ system in the x space direction. 

We consider a structured mesh of size Ax ; -, defined by the cells It = [*j_i/2)*i+i/2)> where we 
have set x ;+1/ / 2 = x/ + Ax/2 ,i G Z, at time t n . As usual, we consider known a piecewise contant 
approximation U h (x,t n ), defined by U h (x,t n ) = U?,x G /, , Vi G Z. 
At initial time t = 0, we impose 



Uf= U Q (x)dx, 

J Xs_ r n 



' x i-l/2 

where Uq is the initial data. This approximation evolves in time, involving a suitable approximate 
Riemann solver. In the HLL approach, the exact Riemann solver solution is substituted by a 
single approximate state (see Figure [2]). Here bi and bR are relevant approximations of Al and 
A/?, respectively. Let us introduce the proposed approximate solution: 

/ vr ( A \ ( ^ L if 7 < bL ' 

Moreover, the search of weak solutions leads to the Rankine-Hugoniot jump conditions 

-b L [W* - W L ) + [F - &{W L )] = 0, (27a) 
-b R [<& R - + [&{W R ) - F] = 0. (27b) 
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L u*f(u l m )b 



U L J(U L ) 




Figure 2: Structure of the approximate HLL Riemann Solver 



These relations provide us with an explicit expression for the intermediate state and flux of the 
numerical scheme 

f = b^m-WW-MHk-*,.) (28b) 

bR - b L 

At each interface x i+ ^/ 2 , we impose the above HLL approximate Riemann solver, assuming the 
CFL like condition (|29l) ensuring that the Riemann solvers do not interact in the case where 
b L ,i+i/2<0andb R j_ l/2 >0: 

A* < b L;l+l/2 b R4 _ l/2 
&x ~ b L j +l / 2 -b R j_ { / 2 

We set fy h (x,t + At), at time t n +At, the superposition of the non-interacting Riemann solutions. 
We define the updated approximation at time t n+l by 



. i r x i+\/i , 

^«+i = / W h (x,t n +At). 
AxJ Xi _ 1/2 

An easy computation gives 

% n+l = % n -^ + i/2-FlU/2), (30) 



where 

r &(% n ) ifo<^ J(+ i /2 

F%. 1/2 W,V&i) = F i+1/2 if b L , i+1/2 < < b R , l+l/2 

( if^m/2<0 

The robustness of the scheme, namely the positivity, the flux limitation, the total mass preserva- 
tion, has been established for the HLL scheme (see (HI for further details). 
Finally, concerning the high order extension, we adopt a van Leer MUSCL technique 11371 , sup- 
plemented by a suitable slope limitation preserving these expected physical properties 0. 
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5.2 An accurate HLLC scheme 



The HLL scheme has proved to be robust, however, its 2D extension fails when approximating 
contact waves. Several works 01 [H introduce a more accurate scheme, the HLLC scheme, based 
on a two state approximation, denoted by U[ and U^. 

First, let us recall the relevant linearization that permits us to define an approximation with 
two intermediate states: on the one hand, the Rankine-Hugoniot conditions (|27l) are considered; 
on the other hand, they are supplemented by the continuity of the Riemann invariants accross the 
contact wave: 



(Ac)l = (AO* = A* > n£ = n£ = n* , (31) 

where /3 X and n are defined by the relation (|24l) . The combination of both the Rankine-Hugoniot 
condition (ITTT) and the relation (1311) standing as the continuity of the Riemann invariants accross 
the contact wave, is sufficient to determine uniquely J3j [H the two approximate states U[ and 
U^, together with their associated fluxes Fl and Fr. The proposed HLLC approximate solution 
can be written as 

r w L iff<&L, 

W*W*t''0 = S f^flff' 02) 

( % if^<f. 

Similar to the derivation of the HLL scheme, we integrate over a cell /, the juxtaposition of the 
non-interacting HLLC Riemann approximate solvers at each interface (projection step), in order 
to obtain the updated quantity 



= _L [ Xi+1/2 W h (x,t n + At). 
AxJ Xi _ 1/2 



This brief description of the HLLC scheme is now completed. It is able to capture exactly the 
contact wave, and satisfies the positivity, the flux limitation, and the total mass preservation. 



6 Numerical Results 
6.1 Central Void 

The first test case is taken from the medical physics literature [2]. We consider only elastic 
scattering, which is modeled by the Henyey-Greenstein kernel. Thus 5 = and 7] n = 0. We 
compare the particle flux \ff(°\x) obtained with the minimum entropy model (labeled Ml) with 
a discrete ordinates solution of the transport equation (labeled SN) with sufficiently many angles 
(128). The method has been described in |[T6ll . 

The test case consists of a one-dimensional geometry with three layers: optically thick, fol- 
lowed by optically thin followed again by optically thick. The layers have an equal depth of 40 
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mm. The scattering and absorption coefficients are o s = 0.5 mm , o a = 0.005 mm for the 
optically thick region, and o s = 0.01 mm , o a = 0.0001 mm -1 for the optically thin region. 
Moreover, g = 0. Figure |3] shows the particle flux as a functon of space. Compared to the 
benchmark solution, the minimum entropy model slightly overestimates but nevertheless quite 
accurately describes the particle flux. In Figure [3] we also show the partice distribution function 
\jf(x,Q.), where Q. = (0,0, /i) in ID. The main difference is that for Mi, the forward-peak of the 
incoming particles reaches further into the medium. 



6.2 Two-dimensional Void-like Layer 

Our second test case, again taken from [O , is a two-dimensional quadratic domain which con- 



tains a void-like layer, shown in gray in Figure |4(a)| Again, we consider only elastic scattering 
modeled by the Henyey-Greenstein kernel. 

We take o s = 0.5 mm 1 and o a = 0.005 mm 1 inside the square, and o s = 0.01 mm 1 
and o a = 0.0001 mm 1 in the void-like ring. In both regions, g = 0. An isotropic source of 
particles is placed on the left boundary. In a 2D contour plot (Figure H]), the fluxes \fr from the 
discrete ordinate method and from the minimum entropy method are virtually indistinguishable. 
The propagation into the medium, as well as the void-like layer are equally well resolved. A 
difference between the models only becomes apparent in a logarithmic plot of a cut through the 
center of the square at y = 50 mm. FigureH] shows the particle flux along this line. The difference 
between both solutions is again of the order of one percent. 



6.3 Electrons on Water Phantom 

As a first test case that includes energy loss, we consider a 10 MeV electron beam impinging onto 
a slab of water. In Figure [5] we compare the results computed with our code to the dose com- 
puted by the state-of-the-art Monte Carlo code PENELOPE [|32l . This code has been extensively 
validated against experimental results. 

To obtain a good fit with the tabulated scattering data, we have fixed our model parameters 
for water as e B = 16.0 eV, Z = 9.40, p el = 0.256 x 10 23 g/cm 3 , p in = 6.21 x 10 23 g/cm 3 . These 
parameters are directly inserted into the model ©> and subsequent derived models issued from 
©. As boundary conditions, we have taken a very narrow Gaussian in energy, and a 8 pulse in 
angle 

Yb = ^oexp(-200(e-e b eam) 2 )5(p- 1), 

and computed the angular moments. PENELOPE was set up in a pseudo-ID setting with a large 
beam size perpendicular to the beam direction. 

In order to compare the different formulations of the models, both depth-dose curves in Fig- 
ure |5] have been normalized to dose maximum one. The penetration depth computed with the M\ 
model agrees very well with the Monte Carlo result. In fact this deviation is within the margin 
of differences between different Monte Carlo codes 11331 . The only major difference occurs near 
the boundary, where the M\ model overestimates the dose. This might be due to the simplified 
physics (possibly neglection of Bremsstrahlung effects) or an oversimplification of the angular 
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(a) Particle flux. 
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(b) Comparison of M\ and benchmark distribution functions. Logarithmic scale. 
Figure 3: Geometry with central void. 
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(c) Transport solution. (d) Minimum entropy model. 



Figure 4: Transport versus minimum entropy for void-like layer. 
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Figure 5: Dose for 10 MeV electron beam on water. 

dependence of i/a in the Mi model. Both possibilities will be investigated further. However, 
we believe that this result can serve as a proof of concept of a PDE based modeling of dose 
computation. 

6.4 CT Data 

In our final test case, we compare our method with Monte Carlo results from PENELOPE using 
real patient CT data showing the hip bone. We took a two-dimensional slice of 6 x 6 cm from the 
three-dimensional CT data. A square region is split into 64x64 squares. In each of the squares, 
the material is described by its Hounsfield grey value (jc,v). The grey values can be translated 
into physical parameters as follows, 

((Six y) \ 
P (X, y)=[ 10Q0 + 1 ) PWater, 

i.e. the densities p e i and p m for water are multiplied by a specified factor. The region shows the 
hip bone and the density varies between 86% and 226% of the value of water. The boundary 
conditions were set up similar to the previous case, with three beams of width 2 cm, each con- 
sisting of 10 MeV electrons, impinging from the centers of three sides of the domain. Contour 
plots of the dose distribution are shown in Figure [7J There, we also show two cuts through the 
dose distribution, looking at the 2D dose distribution, the contour lines agree very well. Note 
that, although we have used 3 x 10 10 particles, there still is significant noise in the Monte Carlo 
results. The two cuts through the beam centers show that also quantitatively the independently 
computed dose distributions agree very well. 
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Figure 6: CT data of hip bone. 
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(c) Cut along x — 3 cm. (d) Cut along y = 3 cm. 

Figure 7: Dose distribution for three beams impinging on hip bone. 
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The computation time for the 3D Monte Carlo dose was 3 x 29 hours for 3 x 10 particles 
on a 3GHz Pentium 4 with 1 GB RAM. In ID, the minimum entropy model took about 1 second, 
in 2D 4 seconds. Thus we expect a computation time of several seconds in a full 3D dose 
computation. 

Again, this result shows that if our model is developed further, it may serve as an alternative 
to existing dose computation methods. 
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